ON THE APPLICABILITY OF CONSTRAINED SYMPLECTIC 
INTEGRATORS IN GENERAL RELATIVITY 



JORG FRAUENDIENER 

Abstract. The purpose of this note is to point out that a naive application of 
symplectic integration schemes for Hamiltonian systems with constraints such as 
SHAKE or RATTLE which preserve holonomic constraints encounters difficulties 
when applied to the numerical treatment of the equations of general relativity. 



It is well known that the equations of General Relativity (GR) can be derived 
from a variational principle and that they can be cast into Hamiltonian form. The 
underlying symplectic structure has been studied as early as the 1940's beginning 
with the work of Bergmann [ , ], Dirac [ , ] and ADM [ ]. The main motivation 
then has been to work out a quantisation scheme for GR. For various reasons, not 
the least of them being the peculiar nature of the symplectic structure of GR, these 
early attempts have not led to any viable theory of quantum gravity. 

On the other hand it has been well established within the numerical mathe- 
matics community [12, 13, It] that the use of so called symplectic integrators i.e., 
numerical ODE solvers which preserve an underlying symplectic structure can 
lead to significant improvements in long-time stability, conservation of first inte- 
grals and accuracy. These methods have been generalised even to Hamiltonian 
systems with constraints. There are two particularly noteworthy methods which 
are called SHAKE [ ] and RATTLE [ ] . They have been developed within the area 
of molecular dynamics but they have since then been used successfully in various 
other applications. However, they only work for holonomic constraints. 

Given the success of these methods it is, therefore, natural to apply symplectic 
numerical methods also to the equations of GR. However, as we will argue in 
this paper, it is not clear (yet) whether there is any advantage to be gained in this 
approach. 

This paper addresses the question of the applicability of symplectic integra- 
tors in GR and it is directed towards both communities, numerical mathematics 
as well as numerical relativists. This necessarily means that we need to review 
both the Hamiltonian framework for GR as well as the essence of symplectic in- 
tegrators. This is reflected in the structure of the paper which consists mostly of 
sections to introduce the necessary background material. In sect. 1 we describe 
finite-dimensional Hamiltonian systems and in sect. 2 we expand this to include 
systems with constraints. Sect. 3 is devoted to a brief exposition of the symplectic 
structure of GR in the special case of spatially compact space-times. In sect. 4 we 
describe the essential properties of symplectic integrators for constrained systems. 
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Finally, in sect. 5 we discuss the consequences of trying to combine these two areas 
of research. 

1. Hamiltonian systems 

Before we come to the symplectic structure of GR let us first look at a classical 
Hamiltonian system with finitely many degrees of freedom such as those occur- 
ring in classical mechanics, molecular dynamics etc. The system is specified by a 
triple (^^, CO, H), where 3^ is a real manifold of even dimension 2n which carries a 
symplectic form co, i.e., a non-degenerate closed 2-form. The pair co) is called 
the phase space of the system. It is the collection of all states which are accessible 
to the system and the symplectic form provides a way to locally sort the degrees 
of freedom into pairs of conjugate variables. 

Since lo is non-degenerate it defines at each point x G ^ an isomorphism be- 
tween the tangent space and the co-tangent space T*^. Thus, any function 
/ E defines a Hamiltonian vector field by the equation 

(1) X|ja; + d/ = 0. 

From this equation and the closure of co follows that the Lie derivative 

(2) Lx^o; = d(X| jo;) + jdo; = 0, 

i.e., the symplectic form is invariant under the flow generated by a Hamiltonian 
vector field. Each member of the flow is a canonical transformation. This is true, 
in particular, for the Hamiltonian vector field X^ generated by the Hamiltonian 
function H : ^ ^ R, the function which specifies the djrnamics of the system; 
the time evolution map <pt generated by X^ which maps an arbitrary initial state 
Xq G >^ to the state at time Hs a canonical transformation. 

Dual to the symplectic form we can introduce a Poisson structure, i.e., Poisson 
brackets {■, ■} on 3^ , by defining for any two functions f,gE 

(3) {f,g}:=a;{Xf,Xg) = Lx^g. 

This turns the algebra of functions on into a Lie algebra with respect to the 
Poisson bracket, the Jacobi identity being a consequence of the closure of a;. It is 
well known that there exist preferred so called canonical coordinates {pi^, q^) on £P 
such that locally the symplectic form is 

CO = dpi( A dq'^ 

or, equivalently, such that these coordinates have canonical commutation relations 

{PuPk} = 0, {ci\q'}=0, {p^,c,'^} = Sf. 

The flow generated by a function H e '^#'°°(^) induces a change in a function / 
which is given by 

In particular, the rate of change in the canonical variables can be used to obtain a 
coordinate expression for the flow 

q' = {H,q'}, Pk={H,pk}- 
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Hamiltonian systems frequently arise from Lagrangian systems by performing 
a Legendre transformation. The most common case is where the Lagrangian sys- 
tem is defined by an action functional 

^ = y ^{q,q)dt 

over a Lagrangian function ^ : TQ ^ ]R on the tangent bundle of a configura- 
tion manifold Q. A Legendre transformation is then used to define a Hamiltonian 
system on the cotangent bundle T*Q of the configuration space. A detailed de- 
scription of these structures can be found e.g., in [1,4,23] 

In many cases the Legendre transformation is well-defined and invertible and 
the Hamiltonian system is valid without any restrictions, i.e., it is unconstrained. 
In some cases, however, when the Lagrangian function is degenerate, the Le- 
gendre transformation is not a local diffeomorphism. This implies that not all 
the possible states in T*Q are available to the Hamiltonian system, i.e., that there 
are constraints which have to be imposed. 

This situation has been analysed in detail by Dirac [ , ] (see also [ ]) who 
developed a theory of Hamiltonian systems with constraints. 

2. Constraints in Hamiltonian systems 

From the geometric point of view a constraint in a phase space a;) is a sub- 
manifold of ^ which comprises the states which are accessible to the system. 
The symplectic form co restricts to a closed 2-form co on'ff. In general, a) will not 
be regular Let 

= {!J e T:,^ : co{U,V) = 0, VV e T/^}. 

At each x G this is a subspace of T^*^ and we assume that the dimension of Gx is 
constant as x varies over Then G = U G^: is a sub-bundle of T'^ (and hence also 
of T^) which defines a distribution in T'^. It is easily seen that this distribution 
is integrable: let X, Y be two sections of G, so that Xjd> = Y jco = 0. Then the 
closure of to implies 

[X,Y] Jd5 = Lx(yja>) -YjLxo) = -Y jd(X jd>) - Y j(X jdd;) = 0. 

Therefore, there exist maximal integral surfaces ^ tangent to G which foliate 
Under certain technical assumptions (for the details see [ ] and references therein) 
the space of leaves = "^[(^ is a differentiable manifold. Furthermore, there ex- 
ists a closed 2-form co' on 3^' which pulls back to co under the canonical projection 
and which is regular. Thus, the pair (^', co') is a phase space on its own. 

This is all that can be said from 'inside 'rf', i.e., without taking into account that 
^ is in fact a sub-manifold of Doing this, one obtains information about how 
the embedding of in ^ affects the structure inside Let us first define 

T/^'-L = {Ue Tx^ ■■ co{U, V) = 0, VV e Tx'^}. 

then, clearly, Gx = Tx'^ia n Tr"^^. Let r be the co-dimension of in then we 
have dimTx'^ = 2n - r and dimT;^'^-^ = r. Furthermore, let / E be 
constant on so that the restriction of df to vanishes. Then at all x E we 
have for any V E Tx^^ 

cVx{Xf,V) = -V{f) = 
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i.e., Xf{x) G Tx'^ia . Let Ca be r independent functions which vanish on near 
X so that '^^ may locally be regarded as the zero-set of these functions. Clearly, 
any function which is locally constant on is functionally dependent on the C/{. 
Hence, the Hamiltonian vector fields '■= evaluated at x generate a r- 
dimensional vector space which, therefore, coincides with Tx^'^. 

The vector fields Xjr for locally constant / need not be tangent to We will 
be interested mostly in two cases: when either none or all of the vector fields are 
tangent to 

In the first case we have Gx = Tx^ (1 Tx^'^ = {0}, so that co is regular and 
3^' = . Then ,Co) is a symplectic sub-manifold of {^S^ ,<xi) i.e., it is a phase 
space in its own right. Note, that Q^g := <jj(Xa, Xg) = {C/i, Cg} is a non-singular 
r X r-matrix when evaluated on . In this case, the constraint functions Ca are 
called second class constraints. 

Since ('^, <i)) is a phase space there exists also a Poisson bracket on corre- 
sponding to d) defined for functions on Denoting the inverse of Qab by Q^^, 
so that QabQ^^ = ^A '-^^ express the Poisson bracket {f,g} between two 
functions / and ^ on in terms of Poisson brackets on ^ as follows. Choose ex- 
tensions of / and g to £P, i.e., functions / and g on ^ which restrict to / and g on 

Then, on "^if the following equation holds: 

(4) {f,g} = {f,g} - {/,C^}Q^«{Cg,g}. 

Here, the left hand side is the Poisson bracket on {'t^,co) and it is defined only on 
while the right hand side is well-defined even an 3^. It makes sense for arbitrary 
fimctions / and g. It is easy to see that it vanishes if / or ^ are taken as constraints. 
Since two extensions of / coincide on they differ by constraints. This shows 
that it is irrelevant which extensions for f or g are used. The expression on the 
right hand side satisfies the defining properties of a Poisson structure so we may 
also regard it as defining a new Poisson bracket {■, jo on which is adapted 
to the existence of the constraint surface. This new Poisson bracket is called Dirac 
bracket [10]. Note, that we can now express the Poisson bracket on in terms of 
Dirac's bracket 

which in turn enables us to discuss the Poisson structure of constrained system in 
terms of quantities on the original phase space. 

The second case of interest is characterised by the fact that all the Hamiltonian 
vector fields Xa corresponding to constraint functions are tangent to . Therefore, 
we have Tx'^^^ C Tx"^^ and Gx ~ Tx'ta^. This implies, that 

{Ca,Cb} = <x;(Xyi,Xg) = Lx^Cg 

which vanishes on It has been useful to introduce the notion of 'weak equality' 
of two functions / and g if and only if they restrict to the same function on 
Thus, 

f^g ^ f-g=}l''CA 

for appropriate functions e "if °°(^). Hence, in the present case we may write 

{CA,Cg}«0. 

In this case, the functions Ca which define the constraint hypersurface are in in- 
volution. They are called /irst class constraints. 
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Since Gx 7^ {0} the restriction of the symplectic form co is degenerate and 
{'rf, co) is a pre-symplectic manifold. Factoring out the leaves of the foliation we 
obtain the reduced phase space {^',lo'), sometimes called the space of the true 
degrees of freedom. 

Let us now consider time evolution. Given a Hamiltonian H G for 
a system with constraints we need to ask for compatibility of the time evolution 
generated by H with the constraints: when the system is started out on '^/f then it 
should remain on '^^ i.e., the Hamiltonian vector field should be tangent to ^ 
or, expressed in terms of Poisson brackets, the weak equality 

(5) {H,Ca}^0 

should hold for all constraints Cj\- Clearly, for the behaviour of the constrained 
system only the restriction H of the Hamiltonian function to is relevant and the 
extensions of J? to ^ (of which H is one) are all a priori equivalent. However, 
we may try to find a compatible extension H for which the Hamiltonian vector 
field is tangent to ^. Writing H = H + A^Cb we find 

« {H,Ca} = {H,Q} + A«{Cb,Q} + {X^CaKb « {H,Q} + A^Qb^. 

This equation tells us that we can find a compatible extension only if Qab is in- 
vertible, i.e., only if the constraints are second class. Only in this case we can 
express the d5mamics of the constrained system entirely in terms of the original 
phase space 

In the case of first class constraints Qab ~ so that either all extensions or none 
satisfy the compatibility condition (5). If it is satisfied then H is constant along the 
Hamiltonian vector fields generated by the constraints Ca ■ Hence, it descends 
to a well-defined function on Furthermore, for its Hamiltonian vector field 
Xh we have 

[Xh,Xa]jcv= -d{{HXA})- 

Since for any weakly vanishing function / ~ one has df = dA^CA + A^dC^^i for 
suitable functions A"^ this implies that for any x and Y e T^^^ 

a;([X^,XH],Y) = co{[Xa,Xh],Y) = A^Y(Cb) = 0. 

Thus, [X^,Xh] e T/^^ so that 

This implies that X^ is projectable onto One can also easily see, that its pro- 
jection is the Hamiltonian vector field for the projected Hamiltonian with respect 
to the symplectic form co'. 

Let us now illustrate the two cases with two examples. 

2.1. Example 1: a particle restricted to a hypersurface. Consider a free particle in 
a Riemannian manifold {Q,gab) whose motion is restricted to a hypersurface S C 
Q. Let Co = F be a function whose zero-set locally defines S. In local coordinates 
q" on Q the action for this situation is given by 



^= J (^'«^«fc('?)4V-AF(^)') dt 
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This leads to the Hamiltonian H = ^g''^PaPb + Requiring that {H, Cq} w 

gives us (using the notation Fa = VoF) 

Ci := p"Fa « 

so we need to include Ci as a constraint. Since {Cq, Ci} = — FaFfc^"* ^ we can 
solve the equations 

{H + A0C0 + AiCi,Q} wOforf = 0,1 
for Aq and Ai and obtain 

-1 ..a ..h ^ 



M =0/ -^0 = 



1 p«p^F„, 



with F^i, = Va^b^- Hence, the final Hamiltonian is 

2m m p'^pa 

It is straightforward to check that its Hamiltonian vector field annihilates both 
constraints. 

2.2. Example 2: relativistic particle. We consider a particle in a Lorentzian space- 
time {Q,gab)- Iri this case the action for the world-line (fir) e Q of the particle is 
given by 



= m J y gabq^cf dr. 

The distinguishing feature of this action is its invariance under reparametrisation, 
T I— > t' = T(t). The conjugate momentum is 

m 

Pa = -Jl gab, 



where we abbreviate £ = \J gab^''^^- Obviously, we obtain the relation 

(6) C{p,q):=g''^paPi-m^ = Q, 

i.e., the momenta cannot attain all possible values. Hence, the states of the system 
are confined to the sub-manifold "^(o C T*Q defined by (6). From this constraint we 
obtain the further relation 

(7) p'dpa = 

which holds on The restriction of a; to has a kernel which we can determine 
as follows. Let X = X''d/dq'' + Y},d/dpj, then we search for non-vanishing X on ^ 
with 

= Xjcj = Y^dq'' - XMpa 
which, in view of (7) implies = and X" = ap" for an arbitrary function a on 
Thus, every vector field in the kernel of a) has the form 

X = ocp"^. 
dq" 

Since the kernel is 1-dimensional the vector fields are proportional to each other 
and their integral curves coincide as sets. It is not difficult to show that these 
vector fields generate exactly the reparametrisation along the integral curves, i.e., 
they generate gauge-transformations. 
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The Hamiltonian vector field of the constraint C is also in the kernel of to 
so that it is tangent to It generates the flow 

The Hamiltonian function can be determined from the Lagrangian in the usual 
way 

tn u 

H = paq" -ml= -jgahq'r -m£ = 0. 



Clearly, this Hamiltonian is compatible with the constraints. In fact, it vanishes on 

which is consistent with the fact that it generates gauge transformations. 

Thus, we have the following picture. The system does not specify individual 
points ipa^q") G as its states but instead one should regard as one state the 
collection of all points which lie on the same integral curve of the gauge vector 
fields X. They must be considered as equivalent because they are related by some 
gauge-transformation. Hence, the states of the system are global entities, an entire 
world-line considered as a point set i.e., without a distinguished parametrisation. 

Since the Hamiltonian vanishes on '^S' it is functionally dependent on the con- 
straint and it also generates a gauge-transformation. So in this sense there is no 
distinguished time evolution in this system which would map from one state to 
another as it is the case in many 'normal' systems. 

If one is interested in the structure of an individual world-line then one can pro- 
ceed by fixing an initial point on the line and then, using the Hamiltonian vector 
field of H, the integral curve through that point can be found. However, the re- 
sult will be a curve together with a special parameter which is determined by the 
choice of the Hamiltonian. The system of a relativistic particle is very similar to 
the situation in GR to which we will now turn. 



3. The symplectic structure of GR 

We now come to a brief introduction to the symplectic structure of GR. We fol- 
low loosely the exposition in [ ]. Other treatments can be found in e.g., [7, 11, 22, 
23]. Let S be a 3-dimensional compact closed manifold^ We consider globally 
hyperbolic space-times of the form M = S x R. We choose a global time-function 
f : S X ]R ^ R and a vector field t" such that the h3^ersurfaces of constant t 
are diffeomorphic to E and such that fd^t = 1. We assume that the hypersurfaces 
Sf are space-like and that the vector field f is future directed and time-like. Let tia 
be the future directed co-normal of the hypersurfaces and denote by ^g„i, resp. g^h 
the space-time metric resp. the metric on . 

We can perform a 3 + 1-decomposition of the geometrical quantities in the usual 
way [ ] by writing f = an" + jH", thereby introducing the lapse function a and 



We concentrate here on the case of spatially closed space-times because we are interested in the 
intrinsic Hamiltonian framework. Issues concerning boundary conditions like in the case of asymptot- 
ically flat space-times or even in the quasi-local regime are somewhat cumbersome to formulate or are 
not even resolved yet [. ]. 
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the shift vector /S". Thus, we can express the 4-geometry in terms of (famiUes of) 
3-dimensional quantities. In this way the Einstein-Hilbert action 

(8) / 'Rj^gd'x 

can be expressed up to boundary terms as the following action 

(9) ^ = / ^{g,g;(^,^)dt 
where the Lagrangian is 

(10) if (g,g; a, ^)^J^a(R + K'^K.^, - K^) ^d^x. 

Here, we have used the scalar curvature R of the metric gab on Sf, the extrinsic 
curvature K^i, and its trace K — Kc'^ of Z,t within the space-time M. Due to the 
relationship 

2KKai, = gab - (L/3^)flfc 

between the extrinsic curvature and the Lie derivative gab '■= 0^tg)ab of the metric 
the Lagrangian is considered as a functional of gab, its time derivative gab as well 
as the lapse and shift. Note, that .if does not contain any time derivatives of a. or 
/S" which indicates that it is singular. In fact, computing the variations of ^ with 
respect to a and /S" yields 

^^^^ C^^ = V^(R-rt«, + x2), Ca^^^^2^Vb{K\-SlK). 

The vanishing of these expressions as required by the Eiiler-Lagrange equations 
yields constraints on the possible configurations. 

In a similar way we compute the momentum conjugate to gab as 

(12) p«i'^^ = ^(r^_Kg«''). 

''gab ^ ' 

Note, that this and the constraint expressions are tensor valued densities of weight 1. 
Finally, we determine the Hamiltonian from the formula 

(13) H{g,p) = IjabP'^dh-^ 
and find (up to boundary terms) 



(14) H{g,p)^fa^ 



+ [-2VaP%] d^x. 



Thus, we have the following situation. As the configuration space Q we take the 
space of Riemannian metrics on Z,. The tangent space TgQ consists of all S5anmetric 
covariant second rank tensor fields 3gab on Z,. The cotangent space T*Q is defined 
as the space of functionally differentiable 1-forms on TgQ, i.e., linear real-valued 
maps which are of the form 



TgQDSgab^ l^p'^Sgab 



where p"^ is a tensor valued density of weight 1. The phase space ^ of general 
relativity (in the context of spatially closed space-times) is the cotangent bundle 
T*Q over the space Q of Riemannian metrics over Z,. Points of ^ are represented 
as pairs {gab'P"^) and tangent vectors to ^ are represented as pairs {Sgab,Sp'^^)- 
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Being a cotangent bundle carries a canonical symplectic form and hence also a 
Poisson structure. 

The symplectic form between two tangent vectors to ^ is defined by 

(15) ^(s,v){{^\?,'hv)\hg,hv)) ^ ^^hv''^hgah-hv"'^hgah<^^ 

and the corresponding Poisson bracket between two functions F and G on ^ is 

(16 |f/G}= / T- 1^ rd'^x. 

The constraints expressions (11) yield functions on ^ by integration over 2L 

Cf = j^fC d^x, Cv = ^ v'Ca d^x, 

where / and v = v" are arbitrary test (vector) fields on S. Using the Poisson 
bracket we can easily see that the constraint functions satisfy the Poisson commu- 
tation relations 





{Cf,Cg} = 




(17) 


{Cf,Cv} = 


-Cv(/)' 




{ Cv/ Cw } ~ 





Therefore, the Poisson brackets among all constraints are again constraints, i.e., 
the constraints are first class. The constraint functions Cj: and Cv generate transfor- 
mations on 't^ which correspond to gauge-transformations, thus mapping a state 
igab'P"^) to 'equivalent' state. The constraints Cy generate 3-dimensional dif- 
feomorphisms within Z. The constraints Cjr, however, generate transformations 
between different hypersurfaces which can be interpreted as the 'evolution' of 
the intrinsic and extrinsic geometry of S within the space-time M along the vector 
field t" = fn". 

The Hamiltonian (14) turns out to be a combination of constraints 
(18) H{g, tt) = Q + C^. 

Hence, it generates gauge-transformations, namely the evolution of S along the 
general evolution vector t" = an" + jS". This implies that we have a similar situa- 
tion here as in the case of the relativistic particle. A particular given state {gabi P"^) 
on is equivalent to states {gabr P"^) which are obtained by such transformations. 
Each equivalence class corresponds to the same single space-time. 

The fact that GR is a completely constrained system is the Hamiltonian way of 
reinstating general covariance of the theory. Any time-evolution in the Hamilton- 
ian sense would map equivalence classes to equivalence classes, i.e., a space-time 
to an entirely different space-time which would not make any sense. Instead the 
Hamiltonian formulation of GR specifies the general covariant geometry of a sin- 
gle space-time eliminating any allusion to a notion of time. 

4. Symplectic integrators 

Let {^,co,H)he a (finite-dimensional) Hamiltonian system possibly with con- 
straints. The flow generated by H maps initial states xq to later states Xf = (^f (xq)- 
The map : ^ — > ^ is a canonical map, the 'time-f map. It is obtained by 
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finding the integral curves of the Hamiltonian vector field of H, i.e., by solving a 
system of ODE when expressed in canonical coordinates. 

There are many methods to solve systems of ODE by numerical means. Some 
of them have the special property that they preserve the structure defining the 
Hamiltonian system. We may regard a numerical method as a map O;, : ^ ^ ^ 
which maps a state x„ to the next state and we call such a method a symplec- 
tic integrator (of order p) if O;, is a canonical transformation for every h which ap- 
proximates the exact Hamiltonian flow for a Hamiltonian function H in the sense 
that 

(19) ^h{x) = cPhix) + i^{hP+') 

for all X G In [ ] it is shown that a symplectic integrator of order p is backward 
stable i.e., that there exists a Hamiltonian function H;, such that Hfj — H = i^{h^'^^) 
and such that O;, is the time-h map of the Hamiltonian vector field corresponding 
to Hf, . This means that a symplectic method can be regarded as the exact time-h 
map for a slightly perturbed Hamiltonian system. 

When constraints are present the symplectic integrators can be generalised to 
numerical methods which preserve the symplectic structure and the constraint hy- 
persurface simultaneously [ , ]. Examples of such methods are the well-known 
algorithms SHAKE [20] and RATTLE [ ] developed within the context of molec- 
ular dynamics. They are implemented schematically as follows. Consider the 
Hamiltonian system (^, co, H) together with constraints C/^ and let xq = {po, qo) 
be a point on the constraint hypersurface We seek a method to compute the 
next point Xj^ = {pi,,qii) after time h on according to the Hamiltonian H. One 
considers the extended Hamiltonian 

H = H + A^Q 

which generates the equations of motion on ^iq 

These have the approximate solutions 

PH = (po - ^^(^o)) - ^'^^^(^o) + ^{h^), 

qH={qo + h^^{xo)^+hA^^-^{xo) + i^{h'). 

However, the multipliers are not yet known. They are determined by requiring 
that the point X;, = (p;„ lies on Thus, one puts 

P = Po-h^{xo), q = qo + h — {xo) 

and notes that 

Cb(x;,) = Cb{£) - ^^(x)^(xo)A^ + ^(x)^(xo)A^ + ^{h^) 

= Cb{£) - hA^ {Cb,Ca} (x) + ff{h^). 
Thus, one can find the multipliers A^ by iteratively solving the linear equation 
(22) CB(x)-/!A^{CB,CA}(i) =0. 
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At each step the are used to update x, thus entering a new iteration until the 
constraints Cb(^) = are satisfied to a desired accuracy at which point one puts 

Due to the special structure of holonomic constraints and their associated 'hid- 
den' constraints the SHAKE and RATTLE algorithms differ in the details of this 
iteration procedure but the general structure of the algorithms is as indicated here. 
The main point about them is that they make the tacit assumption that the matrix 
Qab = {Ca/Cb} is inveriible at every x. This implies that these algorithms work 
only for second class constraints. In fact, the above calculation is nothing but a vari- 
ant of the calculation to find an extension oi H\t^ whose Hamiltonian vector field 
is tangent to 'if. 

5. Conclusion 

We have seen in sect. 3 that GR is a fully constrained theory with first class 
constraints. All the Hamiltonians (14) are combinations of constraints generating 
gauge-transformations. So, strictly speaking, there is no time-evolution. However, 
within computational gravity one uses numerical methods to compute the geome- 
try and hence the physics of one particular space-time. In terms of the Hamiltonian 
framework this can be imderstood as follows. 

Fix initial data, i.e., a point {gab, p"'') on the constraint surface and specify a 
particular Hamiltonian by fixing lapse function and shift vector. This Hamiltonian 
generates a gauge flow which maps the initial point to points which correspond 
to hjrpersurfaces at a 'later' coordinate time. This 'evolution' is clearly symplectic 
and it preserves the constraints. Hence, one can try to use symplectic integrators 
for the task of determining the geometry of the space-time in a particular gauge. 

Suppose that we have arranged a spatial discretisation of the infinite dimen- 
sional system which results in a finite dimensional Hamiltonian system. This 
means that the discretisation results in a system of ODE which is Hamiltonian 
with respect to the discretised symplectic form and which preserves the discre- 
tised constraints. This can be achieved by an appropriate discretization of the ac- 
tion and then performing a Legendre transformation". Let A be a parameter which 
measures the discretisation error. The discretisation should be consistent with the 
continuous system in the sense that we recover the latter from the former in the 
limit A ^ 0. 

Following the implementation of a symplectic integrator we determine the equa- 
tions of motion from an extended Hamiltonian H + \^Cj{. Note, that the index A 
ranges over four times the number of degrees of freedom used in the discretisa- 
tion. As demonstrated in sect. 4 the method relies on the invertibility of the matrix 
Qab = {Ca/Cb}. 

Now two things may happen. Either the Poisson brackets of discretised con- 
straints vanish on the constraint surface i.e., they are also first class with respect to 
the discretised symplectic structure. Then the matrix Qab is not invertible and the 
symplectic integrator algorithm fails. 



It is an interesting and open question as to how much structure of the continuous Hamiltonian 
system can be carried over to the discrete system. 
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The other possibility is that the Poisson brackets of the discretised constraints 
do not vanish which means that the matrix Q^g could be invertible so that multi- 
pliers could be found. However, consistency requires that in the limit of van- 
ishing A one recovers the continuous system from the discrete one. And this in 
turn implies that in that limit the conditioning of the matrix Q^g will become in- 
creasingly bad so that the linear equation (22) cannot be reliably solved anymore. 
Therefore, the continuum limit A — > will result in increasingly inaccurate dis- 
crete approximations to the real solution in contrast to expectations. 

These consequences are observed in numerical implementations of the Einstein 
equations which make use of symplectic integration techniques [ ]. 

The question of how to treat Hamiltonian systems with first class constraints 
numerically appears to be an open issue within the theory of symplectic integra- 
tors. At the moment there is no straightforward remedy to these shortcomings. 
One possibility to circumvent the consequences could be to change the system. 
Recall that we have chosen a Hamiltonian by fixing lapse function and shift vector 
arbitrarily but independently of the evolution. One way to proceed might be to 
couple the choice of these gauge functions to the Hamiltonian system. This could 
break the general covariance in such a way that the resulting system has only sec- 
ond class constraints. However, exactly how to proceed remains largely unclear 
(see [19] for a recent approach). 

Another issue of relevance here is the relationship between holonomic con- 
straints with their hidden constraints on the one hand and the first class / second 
class classification of constraints. Is it possible to find gauge conditions i.e., a 3 + 1 
split and a choice of spatial coordinates, which give second class constraints and 
can be regarded as holonomic constraints in an appropriate generalised sense? 
These issues need further clarifications. 
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